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This  report  summarizes  work  done  between  January  1,  2005  and  December  31,  2006.  The 
original  award  was  made  to  the  University  of  Michigan  (UM)  and  the  duration  was  three 
years,  ending  December  31,  2007.  Dr.  Christlieb  was  the  original  PI,  but  when  he  moved 
to  Michigan  State  University  (MSU)  in  August  2006,  the  contracting  offices  at  AFOSR  and 
the  universities  determined  that  several  modifications  were  needed.  As  a  result.  Dr.  Krasny 
became  the  PI  of  the  award,  and  since  his  portion  of  the  funding  ended  in  December  2006, 
the  award  duration  was  changed  to  two  years,  ending  December  31,  2006.  The  change  in 
the  expiration  date  requires  that  a  final  report  be  written  for  the  UM  award  and  the  present 
document  serves  that  purpose.  The  remaining  funds  were  transferred  to  a  new  one-year 
award  at  MSU  with  Dr.  Christlieb  as  PI,  ending  December  31,  2007.  Hence  the  present  final 
report  corresponds  to  the  first-  two  years  of  the  original  three-year  award  period.  Another 
final  report  will  be  written  for  the  MSU  award,  ending  December  31,  2007. 

1  Objectives 

The  project  aims  to  develop  a  new  grid-free  approach  for  plasma  simulations.  The  specific 
tasks  are  to:  (1)  develop  a  grid-free  field  solver,  fluid  model,  and  kinetic  model,  (2)  evaluate 
these  tools  in  comparison  with  traditional  mesh-based  methods,  and  (3)  demonstrate  the 
capability  of  the  grid-free  approach  in  an  application  of  USAF  interest. 

Most  plasma  simulations  are  currently  performed  using  mesh-based  methods  or  the  hybrid 
particle-in-cell  (PIC)  method  [1,2].  The  alternative  grid-free  approach  developed  here  lias 
the  potential  to  vastly  improve  the  accuracy  and  efficiency  of  these  simulations.  To  achieve 
this  objective  the  investigators  are  applying  several  numerical  techniques  used  in  the  study 
of  vortex  dynamics  in  fluids  [3],  These  techniques  include  (1)  discretization  of  the  flow 
map  by  Lagrangian  particles,  (2)  application  of  a  Cartesian  multipole  treecode  to  accelerate 
the  evaluation  of  electrostatic  forces  induced  by  a  set  of  charged  particles,  (3)  regularizing 
the  Coulomb  potential,  (4)  adaptive  particle  insertion  to  resolve  small  scale  features  in  the 
particle  distribution,  and  (5)  application  of  boundary  integral  techniques  to  impose  boundary 
conditions  on  solid  boundaries. 

'Department  of  Mathematics  -  University  of  Michigan 
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.  2  Status  of  Effort 


As  mentioned  above,  the  present  report  represents  the  first  two  years  of  the  original  award 
period,  up  to  December  31,  2006.  During  that  period,  the  investigators  made  substantial 
progress  in  developing  the  Boundary  Integral/Treecode  (BIT)  approach.  In  particular,  the 
research  effort,  focused  on  several  topics:  (1)  collisionless  plasmas,  (2)  two  stream  instability, 
(3}  crystal  formation  in  trapped  particle  systems,  (4)  collisional  plasmas,  and  (5)  numerical 
heating  in  PIV  vs.  BIT.  The  key  results  are  described  below. 

3  Accomplishments/New  Findings 

The  following  sections  highlight  several  computations  showing  the  feasibility  and  excellent 
performance  of  the  grid- free  approach.  The  method  holds  promise  for  a  variety  of  challenging 
AF  and  civilian  applications  in  plasma  science  and  technology, 

3.1  Collisionless  Plasmas 

The  Vlasov- Poisson  system  describes  a  collisionless  plasma, 

dtfj  +  Vj  •  Vx/j  +  Vv/;  =  0,  (1) 

V^(x.f)  =  (2) 

p(x,t)  =  ^2qjffj(x,v,t)dv,  F ■  = ,  (3) 

j  J 

where  /j(x,  v  J.)  is  the  probability  density  for  a  particle  of  species  j  to  be  located  at  position 
x  with  velocity  v  at  time  t,  (I>  is  the  electrostatic  potential,  q}  is  the  charge  on  species  j,  nij 
is  the  mass  of  a  particle  of  species  j,  p  is  the  total  charge  density,  is  the  permittivity  of 
free  space,  and  F,  is  the  force  acting  on  species  j. 

One  well-accepted  approach  is  to  discretizing  equation  (1)  uses  a  Lagrangian  framework 
in  which  /(x.v.  f)  is  represented  as  a  collection  of  macro-particles  evolving  in  phase  space 
by  Newton’s  equations, 


J7lj 


For  the  electrostatic  problem  in  this  case,  the  inter-particle  force  has  the  form, 


N 

f,  = 

k= 0 
M» 


(4) 

(5) 


Particle-in-cell  (PIC)  is  a  well  known  method  for  evaluating  F,  in  0(N  log  N)  operations 
using  a  fixed  regular  mesh  and  a  fast  Poisson  solver  [2],  However,  PIC  suffers  from  a  host  of 
mesh-based  effects  that  can  significantly  degrade  the  quality  of  the  solution.  These  include 
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<  numerical  heating,  difficulty  with  complex  domains,  and  problems  resolving  internal  and 
boundary  layers  [1]. 

The  alternative  approach  developed  here  couples  a  Lagrangian  description  of  the  Vlasov 
equation  (1)  with  an  efficient  grid-free  field  solver  for  the  Poisson  equation  (2).  The  grid-free 
field  solver  employs  a  boundary  integral  method  and  a  treecode  algorithm  to  evaluate  the 
field  in  time  0(N  log  N)»  The  basic  procedure  is  to  cast  the  Poisson  equation  (2)  as  an 
integral  equation,  using  Green’s  theorem,  and  to  then  discretize  the  resulting  volume  and 
surface  integrals  in  such  a  way  as  to  make  the  system  equivalent  to  a  sum  over  N  charged 
particles.  The  flow  map  is  used  to  change  variables  and  rewrite  the  volume  integral  as  an 
integral  over  a  set  of  Lagrangian  parameters.  Accuracy  is  maintained  by  inserting  new  phase 
space  particles  using  interpolation  with  respect  to  the  Lagrangian  variables.  Initial  timing 
results  show  that  for  high  accuracy  (less  than  a  1%  error  in  FT),  the  BIT  approach  is  superior 
to  a  mesh-based  PIC  method. 

We  applied  the  BIT  approach  to  several  plasma  problems  that  exhibit  complex  dynamics. 
Two  examples  are  the  two  stream  instability  in  neutral  plasmas  and  crystal  formation  of  a 
single  charged  species  confined  in  a  Penning-Malmberg  trap.  These  examples  are  described 
next. 

3,2  Two  Stream  Instability 

The  two  stream  instability  resembles  a  Kelvin- Helmholtz  instability  in  phase  space.  The 
basic  premise  is  that  given  a  neutral  plasma  with  two  electron  streams  flowing  in  opposite 
directions,  in  the  presence  of  a  slight  perturbation,  the  streams  interact  in  a  nonlinear  way 
so  as  to  roll  up  in  phase  space.  The  resulting  vortices  represent  the  trapping  of  electrons  in 
a  phase  space  volume. 

Figure  1  shows  the  evolution  of  the  phase  space  instability  using  the  BIT  approach.  Pe¬ 
riodic  boundary  conditions  in  space  are  imposed  in  this  test  case.  It,  is  worth  noting  that  the 
method  maintains  symmetry  on  very  long  time-scales  whereas  mesh-based  methods  typically 
lose  symmetry  quickly  due  to  interpolation  errors  [2].  Furthermore,  the  BIT  solution  differs 
substantially  from  the  behavior  seen  in  PIC  simulations,  e.g.  the  phase  space  distribution 
does  not  approach  a  Gaussian  as  predicted  and  accepted  by  the  PIC  community.  We  expect 
BIT  to  display  different  behavior  from  PIC  because  we  are  able  to  maintain  much  higher 
accuracy.  Currently  we  are  studying  the  behavior  of  solutions  wit  h  initial  finite  temperature 
distributions.  The  work  on  the  warm  two  stream  Instability  has  involved  the  development  of 
volumetric  point  insertion  using  a  bilinear  interpolant,  which  extends  our  previous  approach 
in  which  particles  were  inserted  along  curves.  These  simulations  were  done  in  collabora¬ 
tion  with  a  talented  undergraduate,  Benjamin  Sonday,  who  is  now  a  graduate  student  at. 
Princeton  University  with  a  DOE  Computational  Science  Graduate  Fellowship.  Details  of 
volumetric  point  insertion  will  appear  in  a  forthcoming  article. 

3-3  Crystal  Formation  in  Trapped  Particle  Systems 

A  Penning-Malmberg  trap  is  a  strongly  magnetized  plasma  in  which  perturbations  in  the 
magnetic  field  due  to  the  flowing  plasma  can  be  neglected  [4].  The  trap  is  a  grounded 
evacuated  metal  cylinder  with  a  large  magnetic  field  in  the  axial  direction,  which  provides 
radial  confinement,  via  the  Lorenz  force,  and  end  cap  potentials  serve  to  contain  a  single 
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Figure  1:  Two  stream  instability.  Plot  of  particle  distribution  in  phase  space,  left:  t  =  0, 
center:  t.  =  2.5,  right:  t  —  3.3.  This  simulation  used  adaptive  point  insertion  to 
maintain  resolution  as  the  particle  beams  roll  up. 


charged  species,  The  system  has  a  large  separation  of  time  scales  and  can  be  reduced  to  the 
following  2D  description,  as  long  as  the  system  is  collisionless. 


dtne  4-  vc  •  Vxnc  =  0 


— V$  x  B 

ImT 


(G) 


where  nc  is  the  electron  density  in  the  trap,  B  is  the  magnetic  field,  <I>  is  the  solution  to  Hit* 
Poisson  equation  (2)  with  $  —  0  on  the  domain  boundary,  and  the  flow  is  in  the  :n/-plane. 
Using  the  grid-free  BIT  method  with  point  insertion,  we  have  studied  crystal  formation 
in  this  system.  Figure  2  shows  a  comparison  between  experiment,  and  simulation.  The 
setup  is  seven  dense  patches  of  electrons  placed  in  a  large  diffuse  background.  The  process 
of  crystal  formation  in  this  system  is  a  dynamic  process  which  depends  on  an  interplay 
between  the  background  and  the  strong  patches.  Numerical  studies  have  shown  that  without 
the  background,  the  system  will  not  crystalize.  Furthermore,  the  system  exhibits  sensitive 
dependence  on  initial  conditions.  The  simulations  show  that  slight  changes  in  the  density 
ratio  and  radius  of  the  initial  dots  greatly  alter  the  critical  time  scales  in  the  system  as 
well  as  the  structure  of  the  crystalline  state.  For  initial  conditions  consistent  with  the 
experiment,  the  simulations  are  in  good  agreement  with  the  experimental  data,  as  shown  in 
Figure  2.  Related  grid- free  results  can  be  found  in  [5  7].  We  plan  to  continue  this  work  by 
implementing  volumetric  point  insertion  in  our  Penning  trap  code. 


3.4  Collisional  Plasmas 


To  account,  for  collisions,  the  right  hand  side  of  equation  (1)  becomes 


(1L 

6 1 


c 


f  [  (f(y')fW  ~  /(v)/(v.))  vadQ  r/v,, 
je.3  Js 


(7) 


where  v  is  the  relative  velocity  between  two  particles,  a  is  the  scattering  cross  section  and 
dQ  is  the  differential  solid  angle  through  which  particles  scatter.  The  two  product  terms  in 
t  he  parentheses  represent  a  source  and  loss  of  particles  which  scatter  into  and  out  of  a  given 
velocity  element  through  collisional  events.  This  convolution  integral  over  velocity  space  is 
expensive  to  compute  and  one  common  approach  to  reducing  the  cost  is  to  use  Monte  Carlo 
techniques  [8]. 

In  this  work  we  have  developed  an  innovative  grid- free  approach  to  Monte  Carlo  methods 
for  gas  dynamics  problems.  Again,  we  model  /  as  a  collection  of  macro  particles  in  phase 
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Figure  2:  Crystal  formation  in  a  Penning- M al mb e rg  trap.  From  left  to  right,  the  panels 
are  in  order  of  increasing  time.  The  top  row  is  experimental  data  taken  at  the 
University  of  California  San  Diego  [4j.  The  second  and  third  rows  are  two  methods 
of  visualizing  the  same  BIT  simulation  results.  Row  two  is  a  visualization  of  the 
actual  Lagrangian  points  and  row  three  uses  a  weighting  t.o  a  mesh  which  is 
similar  to  the  experimental  data  analysis  procedure. 


space.  The  evolution  of  the  new  equation  (1)  is  carried  out  using  operator  splitting  in  time, 
i.e.  we  evolve  the  macro-particles  which  describe  the  transport  of  /  through  phase  space 
and  then  correct  the  resulting  distribution  functions  for  collisional  events  using  a  Monte 
Carlo  method.  It.  lias  been  shown  that  this  method  converges  to  the  solution  of  the  modified 
equation  (1)  as  the  particle  number  goes  to  infinity  [9],  Novel  aspects  in  this  work  are 
that  particles  are  clustered  adaptively  using  octrees  at  each  time  step  and  that  macro  flow 
quantities  are  tracked  and  updated  in  the  nodes  of  the  tree  using  grid-free  interpolation. 

The  method  has  been  applied  to  several  benchmark  problems  such  as  Couette  flow, 
thermal  Couette  flow,  and  flow  past  a  lifting  body.  In  addition,  we  have  applied  the  method 
to  study  the  formation  of  a  Bose-Einstcin  condensate.  This  work  is  done  in  collaboration 
with  a  PhD  student,  Spencer  Olson.  An  article  is  being  prepared  for  submission. 

3.5  Numerical  Heating  in  PIC  vs.  BIT 

The  focus  in  this  work  is  on  the  simulation  of  high  density  plasmas,  as  in  laser-plasma 
interactions.  In  this  case,  PIC  simulations  exhibit  numerical  heating  when  using  compu¬ 
tationally  feasible  mesh  sizes.  We  show’  how  this  artifact  can  he  reduced  with  a  different 
choice  of  time-integration  scheme.  The  main  conclusions  of  the  study  may  be  summarized 
as  follows.  In  both  BIT  and  PIC,  numerical  heating  can  be  reduced  by  employing  implicit 
Euler  (IE)  time  integration.  In  the  PIC  scheme  using  standard  leap-frog  (LF)  time  stepping, 
the  solution  appears  converged  with  respect  to  the  time  step,  but  on  switching  to  IE,  furt  her 
substantial  improvements  are  achieved.  BIT  and  PIC  perform  equally  well  using  IE  or  LF. 
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when  PIC  uses  a  very  fine  mesh,  but  BIT  does  even  better  as  the  time  step  is  reduced. 
Since  BIT  does  not  employ  a  mesh,  it  does  better  than  PIC  when  PIC  uses  a  coarse  mesh. 
In  a  situation  with  too  few  particles  in  the  Debye  sphere,  PIC-LF  code  can  do  better  than 
BIT-LF  due  to  the  spatial  filtering  induced  by  the  mesh,  but  BIT-IE  is  more  accurate  (with 
lower  numerical  heating)  than  PIC-IE.  We  reiterate  that  the  importance  of  this  study  is  the 
implication  for  dense  laser-plasma  problems.  This  work  was  carried  out  by  Dr.  Christlieb  as 
part  of  an  NRC-AFOSR  summer  faculty  fellowship  at  AFRL  Edwards  in  collaboration  with 
Dr.  Jean-Luc  Cambier. 
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doc  support. 

•  AFOSR,  FA9550-06- 1-0529.  09/01/06-05/31/07.  $50,000,  “A  Grid-Free  Particle  Method 
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Invited  Paper 


Abstract — A  common  approach  to  mod  ding  kinetic  problems  in 
plasma  physics  is  to  represent  the  plasma  as  a  set  of  Lagrangian 
macro- particles  which  interact  through  long-range  forces.  In  the 
well-known  particle-in-cdl  (PIC)  method,  the  particle  charges 
are  interpolated  to  a  mesh  and  the  Helds  are  nhtained  using 
a  fast  Poisson  solver.  The  advantage  of  this  approach  is  that 
the  electrostatic  forces  can  be  evaluated  in  time  O(JVlogJV), 
where  N  is  the  number  of  macro-particles*  hut  the  scheme  has 
difficulty  resolving  steep  gradients  and  handling  nonconforming 
domains  unless  a  sufficiently  fine  mesh  is  used.  The  current  work 
describes  a  grid-free  alternative,  the  boundary  integral/!  reecode 
(BIT)  method.  Using  Green’s  theorem,  we  express  the  solution  to 
Poisson’s  equation  as  the  sum  of  a  volume  integral  and  a  boundary 
integral  which  are  computed  using  particle  discretizations.  The 
treecode  replaces  particle-particle  interactions  by  particle-cluster 
interactions  which  are  evaluated  by  Taylor  expansions.  In  addi¬ 
tion,  the  Green's  function  is  regularized  and  adaptive  particle 
insertion  is  implemented  to  maintain  resolution,  Like  PIC,  the 
operation  count  is  0(N  log  N),  but  BIT  avoids  using  a  regular 
grid,  so  it  can  potentially  resolve  steep  gradients  and  handle 
complex  domains  more  efficiently.  We  applied  BIT  to  several 
hounded  plasma  problems  including  a  one-dimensional  fl-D) 
sheath  in  direct  current  (dc)  discharges,  1-D  virtual  cathode,  cold 
two-stream  instability,  two-dimensional  (2-1) )  planar  and  cylin¬ 
drical  ion  optics,  and  particle  dynamics  in  a  Penning-Malmherg 
trap.  Some  comparisons  of  BIT  and  PIC  ivcre  performed.  These 
results  and  ongoing  work  will  be  reviewed. 

Index  Terms — Boundary  integral  method.  Coulomb  potential, 
grid-free,  ion  optics,  multipole  expansion,  partide-in-cell  method, 
Penning-EVI a  1  m he rg  trap,  Poisson  solver,  sheath  formation, 
treecode  algorithm,  two-stream  instability,  virtual  cathode. 


L  Introduction 

MANY  important  problems  in  plasma  physics  involve 
hounded  domains  with  complex  geometry,  e.g.,  space¬ 
craft  thuster  plume  interactions  [1],  plasma  sensors,  and 
semiconductor  fabrication  systems  [2],  [3].  Plasmas  also  arise 
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in  unbounded  domains  in  the  space  environment  and  give  rise 
to  internal  layers  (e.g.,  solar  flares).  This  paper  reports  on  the 
development  of  a  novel  grid-free  approach  for  simulating  these 
type  of  problems. 

Nonequilibrium  plasmas  can  be  described  by  a  variety  of 
models  depending  on  the  physical  conditions.  Complete  infor¬ 
mation  is  obtained  by  solving  the  Boltzmann  equation  for  the 
distribution  function  v,}  of  each  plasma  species  i 


Ofi 

Of 


+  vi  •  Vx./i  +  —  ■  VJ, 

m; 


(1) 


where  x,  is  position,  vp  is  velocity,  F,  is  the  applied/self  force, 
mp  is  mass,  and  {0fi/0t}c  is  a  convolution  integral  describing 
particle-particle  collisions  [4].  Because  (1)  is  six-dimensional 
plus  Lime,  il  is  not  computationally  feasible  lo  solve  ii  for  each 
species.  However,  lo  make  accurate  predictions  about  ihe  be¬ 
havior  of  the  system*  a  kinetic  description  is  often  necessary  for 
at  least  one  of  the  species.  For  example,  the  tail  of  the  electron 
distribution  function  (EDF)  plays  a  dominant  role  in  the  ioniza¬ 
tion  process.  Often  the  energetic  tail  is  far  from  a  Maxwellian 
distribution,  hence,  a  kinetic  model  is  necessary  to  accurately 
describe  the  EDF  and  thereby  correctly  describe  the  ionization 
process  [5],  [6]. 

Many  options  are  available  for  the  numerical  solution  of  the 
Boltzmann  equation  (1)  [7]— [  1 5 J.  In  the  case  of  collisionless 
plasma,  (I)  reduces  to  the  Vlasov-Poisson  (VP)  system.  One 
well-accepted  approach  is  solve  the  VP  system  using  a  La- 
grangian  framework  (LF),  In  an  LF,  /(fx.v)  is  represented 
as  a  collection  of  macro-panicles  in  phase  space  [15],  The 
particles  are  given  an  initial  distribution  by  sampling  phase 
space  at  t  =  0  and  the  evolution  of  the  system  is  obtained  by 
solving  Newton’s  equations  for  each  particle 


In  this  approach,  collisions  can  be  accounted  for  by  coupling 
this  method  to  a  Monte  Carlo  technique  [15]. 

The  electrostatic  force  in  three  dimensions  is  given  hy 


N 


i-n 


(3) 


for  j  —  1 .....  jV.  Direct  summation  requires  ()(N2)  opera¬ 
tions,  hut  Fj  can  he  evaluated  in  0(N  log  N)  operations  by  the 
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inter- particle  distance  (x/Ax) 


Fig,  I  Coulomb  force  as  a  function  of  distance  between  two  charged  particles 
within  a  mesh  cell;  one,  iwo,  and  three  dimensions  compared  with  linear 
weighting  PIC  1 151. 

partide-in-ceil  (  PIC)  method.  PIC  involves  four  stops:  panicles 
are  first  interpolated  to  a  mesh  to  give  a  density,  a  fast  Poisson 
solver  is  used  to  compute  the  fields,  the  force  is  interpolated 
back  to  the  particle  locations,  and  finally  the  new  x,  and  v,  are 
computed  [7].  [8]. 

PIC  has  optimal  performance  when  the  underlying  mesh  is 
uniform.  This  is  because  on  a  nonuni  form  mesh,  standard  finite 
difference  stencils  become  lower  order  accurate  and  more  elabo¬ 
rate  stencils  are  required  to  restore  higher  order  accuracy.  Other 
mesh-based  effects  include  difficulty  in  resolving  internal  layers 
1 16)  and  an  incorrect  description  of  imerpartiele  forces  within  a 
mesh  cell  [  15].  Moreover,  not  resolving  imerpartiele  forces  can 
have  a  significant  impact  on  systems  where  Coulomb  collisions 
are  important,  as  in  high-density  plasmas  [  1 7).  To  illustrate  this. 
Fig.  1  plots  the  exact  Coulomb  force  as  a  function  of  distance  be¬ 
tween  two  charged  panicles  within  a  mesh  cell  in  one,  two,  and 
three  dimensions,  compared  to  the  result  from  a  linear  weighting 
PIC  scheme.  The  PIC  results  differ  considerably  from  the  exact 
force.  This  can  lead  to  spurious  plasma  healing,  but  the  opposite 
occurs  in  a  damped  system,  i.e.,  the  plasma  can  be  erroneously 
cooled  [7].  To  see  why  this  is  a  problem,  consider  a  case  in  which 
scattering  due  to  small  angle  electrostatic  collisions  is  impor¬ 
tant-  In  a  standard  PIC  scheme,  the  imerpartiele  force  fails  to 
zero  inside  a  mesh  cell  (Fig.  I ),  which  is  unphysical.  The  small 
angle  scatter  is  then  computed  incorrectly,  leading  to  unphy steal 
cooling  of  some  particles. 

Much  effort  has  gone  into  overcoming  these  difficulties.  Fi¬ 
nite  element  PIC  addresses  the  issues  of  complex  domain  ge¬ 
ometry  and  boundary  layers,  but  there  arc  still  problems  with  in¬ 
ternal  layers  and  resolving  interparticle  forces  within  a  mesh  cell 
[18],  [19].  The  partide-particle/partide-mesh  (F3M)  method 
was  developed  to  address  the  issue  of  imerpartiele  forces  [8], 
but  the  cost  approaches  0(jV2),  Another  method  for  handling 
short-range  Coulomb  forces  sums  the  collisions  into  a  single 
large  angle  event  [20],  [211,  but  resolving  sharp  gradients  re¬ 
mains  an  issue.  The  current  state  of  the  art  for  plasma  simu¬ 
lations  combines  PIC  with  adaptive  mesh  refinement  (AMR) 
[22]— [25 ],  but  despite  the  success  of  AMR-PIC,  some  limita¬ 
tions  persist.  For  time-dependent  simulations  in  complex  ge¬ 
ometry,  AMR  requires  elaborate  meshing  tools  to  generate  the 


locally  refined  meshes.  Moreover,  in  regions  where  coarse  and 
fine  meshes  intersect*  the  change  in  resolution  can  lead  to  spu¬ 
rious  features.  For  hyperbolic  transport  equations,  this  is  a  local 
defect  that  can  he  corrected,  hut  for  elliptic  field  equations  the 
difficulty  is  nonlocal  and  a  small  error  in  one  location  caji  affect 
the  entire  solution.  Solutions  to  this  problem  involve  filtering 
[25],  but  this  can  introduce  other  ant  facts. 

Here,  we  describe  an  alternative  grid-free  approach  for 
plasma  simulations  in  complex  domains,  the  boundary  in- 
teg  ra  l/l  reecode  (BIT)  method.  Treecodes  were  introduced 
as  an  efficient  method  for  computing  gravitational  forces  in 
systems  of  point  masses  [26 f  The  key  idea  is  to  replace  the 
part  tele-particle  interactions  by  particle-cluster  interactions. 
In  the  simplest  case  of  a  monopole  approximation,  a  cluster 
of  panicles  is  replaced  by  a  single  panicle  carrying  the  total 
mass  of  the  cluster.  This  reduces  the  work  needed  to  compute 
long-range  forces  from  G(N2)  to  ()(N\ogN).  To  gain  effi¬ 
ciency,  the  fast  multi  pole  method  (FMM)  uses  higher  order 
mul  ti  pole  expansions  and  cluster-cluster  interactions  [27], 
[28],  while  other  approaches  use  variable  order  expansions  and 
locally  adapted  clusters  [291. 

Treecodes  and  the  FMM  have  been  used  extensively  in  as¬ 
trophysics  [30],  fluid  dynamics  [31],  and  molecular  dynamics 
[32],  Treecode  simulations  have  also  been  conducted  fur  dense 
plasmas  in  simple  domains  w  ith  free-space  or  periodic  boundary 
conditions  [17],  [33],  [34],  However,  realistic  plasma  systems 
generally  require  Dirichlet,  Neumann,  or  mixed  boundary  con¬ 
ditions  on  a  domain  with  possibly  complex  geometry.  In  this 
case,  a  boundary  integral  method  can  be  used  to  correct  the 
free-space  potential  and  account  for  the  presence  of  a  boundary 
while  retaining  the  essential  grid- free  nature  of  the  scheme  [351, 
[36].  Boundary  integral  methods  use  Green's  formula  and  recast 
the  original  partial  differential  equation,  say  Poisson's  equation, 
as  an  integral  equation  [38],  which  reduces  the  dimensionality 
of  the  problem  by  one. 

The  grid- free  BIT  approach  is  an  attractive  alternative  to 
mesh-based  PIC  methods.  One  advantage  of  a  BIT  approach  is 
that  it  eliminates  errors  introduced  by  interpolation  to  the  mesh. 
Moreover,  BIT  can  naturally  handle  systems  with  nonuni  form 
particle  density  and  complex  geometry.  We  have  applied  BIT 
to  simulate  one -dimensional  (1-D)  sheath  formation  in  direct 
current  (dc)  discharges,  1-D  formation  of  a  virtual  cathode, 
two-stream  instability,  two-dimensional  (2-D)  planar  and 
cylindrical  ion  optics,  and  charged  particle  dynamics  in  a 
Pcnning-Malmherg  trap.  Here,  we  present  preliminary  timing 
comparisons  between  BIT  and  mulligrid  PIC  for  2-D  planar 
geometry.  In  the  following  sections,  we  outline  the  BIT  method, 
discuss  its  advantages  and  disadvantages,  present  a  summary 
of  these  applications,  and  conclude  with  a  discussion  of  future 
prospects. 

II.  BIT 

The  BIT  scheme  combines  a  grid- free  field  solver  with  a 
boundary  integral  method,  and  we  may  view  the  approach  from 
either  of  two  perspectives.  First,  BIT  is  complementary  to  PIC 
in  Lhe  sense  that  the  mesh-based  field  solver  in  PIC  is  replaced 
by  a  grid-free  field  solver  in  BIT.  The  second  perspective  is  to 
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view  the  Lagrangian  macro -part ides  as  marker  points  on  phase 
space  curves-  As  the  curves  evolve,  the  marker  points  spread  out, 
and  new  points  are  inserted  along  the  curves  to  maintain  resolu¬ 
tion.  Point  insertion  is  performed  by  interpolation  with  respect 
to  a  Lagrangian  parameter  along  phase  space  curves  [29],  [37 J. 
To  accommodate  point  insertion,  the  present  formulation  of  the 
field  equation  is  more  general  than  in  our  previous  work  [35], 
[36],  In  the  following  sections,  wc  derive  a  Green's  function 
method  for  Poisson’s  equation,  convert  the  volume  integral  to  a 
sum  over  Lagrangian  macro-particles,  show  how  multipole  ex¬ 
pansions  can  be  used  to  evaluate  the  sum  efficiently,  and  demon- 
si  rate  how  to  handle  boundary  conditions- 


respectively.  Letting  h  —  be  the  solution  of  Poisson's  equa¬ 
tion  (4)  and  ;/  =  G  the  free-spaee  Green’s  function,  we  express 
the  potential  as 


Z\ 

*(#)  =  [  -G{x\y)^dx 
J 


+  \^(x)^G(x\y)-G{r,\u)~^{X) 


(7) 


and 


A,  Potential  Function 
Poisson’s  equation  is 


V 2  <1*  ( x )  =  a$(z)  +  ftdn^[z)  =  7(z)  (4) 

where  4>(x)  is  the  electrostatic  potential  at  xeO  \  Bi L  p(x)  is 
the  charge  density,  is  the  permittivity  of  free  space,  0n  is  the 
normal  derivative,  z  €  and  the  constants  ft,  h  determine  the 
type  of  boundary  condition  (Dirichlet,  Neumann,  or  mixed).  To 
solve  the  Poisson  equation,  we  employ  Green’s  formula 


and 


*(y)  =  /  -G<x|y)^rfx 

n 

+  /(*<x)V„G(x|y)  -G(x|y)V»*(x))  ■  n,lSx  (8) 
on 

where  {y>y}  €  U\  BQ.  Note  that  <I>(y)  =  4*p(y)  +  4J//(y), 
where  the  volume  integral,  is  the  particular  solution  of 

(4),  and  the  boundary  integral,  is  the  homogeneous  so¬ 

lution.  For  later  reference,  we  record  the  particular  solutions, 

Z  L 

$p{y)  ~  -  f  G(x\y)?^dx  (9) 

J  fo 

-0 

and 

$p(y)  =  -  /c(x|y)— -tlx  (10) 

J 

o 

and  the  homogeneous  solutions 


I {hV*g  -  r,V2h)dx  =  i{hVg  -  yVh)  ■  n  dS 
n  on 

in  one  and  two  dimensions,  respectively.  Note  that  the  bold  dx 
implies  a  2-D  integral  over  the  coordinate  x  £  Q.  We  also 
use  the  free-spaee  Green’s  function  for  the  Laplace  operator, 
G(x|y),  which  satisfies 


VyG(x|y)  =  ©5  (||x  -  y||j) ,  [0nG]y  =  G,  [G)y  =  0 


<M.v)  = 


®(x)-jjG(x\y)  -  G{x\y)-£-9{x) 


(ID 


and 


where  ||x  -  y||2  -  “  Vi)2*  [  ]y  denotes  a  jump  at 

y,  and  C  depends  on  the  dimension.  The  free-spaee  Green's 
functions  in  one  and  two  dimensions  are 


<My)  =  j  (<l>(x)V„G(x|y)  -  G(x|y)V*$(x))  ■  ndSx. 
on 

(12) 

Depending  on  the  boundary  conditions,  4>w(y)  can  be  repre¬ 
sented  as  either  a  single  layer  or  a  double  layer  potential.  This 
will  be  discussed  below. 

B ,  Volume  In tegra l/Pa rtic ti la r  Sol u lion 

Recall  that  the  charge  density  in  general  is  given  by 


G(x\y) 


~  V\ 

2 


and 


G{x|y)  =  - 


tB(llx-ylla) 

2  77 


(5) 


P(x) 


v}f/v 


(13) 


where  d  is  the  number  of  velocity  dimensions  and  m  is  the 
number  of  species.  We  consider  two  approaches  for  obtaining 
W  the  potential  <!*/>, 
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*  Particle  Approximation:  In  our  first  approach,  we  represent 
the  phase  space  distribution  /*  as  a  set  of  della  functions 


Ni 

/;  ( / .  X .  V )  =  Vj) 

i- i 

where  w}  is  die  weight  at  {x, .  v, },  Substituting  this  into  the 
equation  for  p  gives 


Fig.  2.  Flow  map  in  phase  space  tc  <  Initial  lime.  T :  Later  time. 


m 

j,(x)  =  ^  ~  xj) 

i=l  ;= 1 

and  the  particular  solution  then  reduces  to 


m  Nt 

My)  =  “  5Z  S  ~^G(xjly)-  (14) 

f-ij-i  fn 

This  is  the  same  as  in  PIC,  where  the  phase  space  particles  x_, 
are  viewed  as  macro -particles. 

integral  Form :  In  our  second  approach,  we  restrict  attention 
to  the  Vlasov  equation  so  that  the  collision  operator  in  the  Boltz¬ 
mann  equation  (1)  is  set  to  zero.  We  make  this  assumption  in 
order  to  focus  on  prohlems  in  which  diffusive  effects  due  to  col¬ 
lisions  are  negligible.  In  the  absence  of  such  diffusive  effects, 
the  solution  of  the  Vlasov  equation  can  develop  fine  scale  fea¬ 
tures  and  these  are  handled  here  by  point  insertion.  Extending 
this  approach  to  colli  sional  plasmas  is  a  topic  for  future  work. 

In  this  approach,  we  employ  ideas  from  Lagrangian  vortex 
methods  in  computational  fluid  dynamics  [31],  The  volume  in¬ 
tegral  over  x  for  ‘bp,  given  above,  is  recast  as  a  volume  integral 
over  coordinates  (x0,  vn),  where 


x  =  x(t,x„v0),  V  =  v(/.x„.v0) 

is  the  flow  map  in  phase  space,  Le.,  (x„.  vrJ)  is  the  initial  con¬ 
dition  (location,  velocity)  that  maps  into  (x,  v)  at  a  later  time  t 
The  initial  condition  (xOJ  vt))  is  a  Lagrangian  parameterization 
of  phase  space  and  this  is  critical  to  the  point  insertion  scheme, 
as  explained  below.  The  volume  integral  over  initial  condition 
for  is  then  discretized  using  the  trapezoidal  rule. 

Before  recasting  the  integral  for  it  helps  to  think  about 
how  we  treat  the  distribution  function  /, .  Consider  a  volume  Ff 
in  phase  space  at  the  initial  time  t  —  t()  and  a  later  time  t  —  T 
as  in  Fig,  2,  Note  that  the  sketch  is  intended  to  represent  one, 
two,  or  three  dimensions  in  space  and  velocity.  Using  the  change 
of  variables  formula,  the  integral  of  /,  over  volume  I>  can  be 
expressed  as  an  integral  over  the  initial  volume 

jj  /;(7\x,  v)dxrtv 

*  Fr 

=  ff  fi(T,x(T,x„,v„),v(T,xn'Vn))J0dx0dv0  (15) 

'  r'>„ 


where  Jn  —  |0{x,  v)/0(xn,  V0)|  is  the  Jacobian  determinant  of 
the  flow  map.  It  is  well  known  that  if  the  collision  operator  is 
identically  zero,  the  fraction  of /r  in  a  phase  space  volume  I\, 
is  preserved  under  the  action  of  external  fields,  and,  therefore, 
in  the  absence  of  absorbing  boundaries,  the  volume  integral  of 
/,  over  phase  space  must  be  constant  for  all  T  >  tti 


f  f  fi(T,x,v)dxdv  =  f  f 
rr  r,„ 


/,(*n.xn,  v„)r/x„r/v(1. 


(16) 


From  (15)  and  (16),  it  follows  that 


fi  (T,  x(T,  x,>,  v0),  v(r.  x(,.  vrt))  ./,>  =  vM)  (17) 

which  expresses  the  conservation  of  charge  in  phase  space  for 
the  Vlasov-Poisson  equation. 

Now  recall  the  particular  solution  in  (10).  Substituting  the 
charge  density  from  (13)  into  thai  expression  gives 


My)  =  -V//i  f  f  G^yYi(r.x.v)rfxrfv. 

-=i  v  fn 


Under  our  change  of  variables,  the  integrand  is 

V°)ly) /.  (T.  x(T.  x0,v„),  v(T.  x„,  v„))  J„ 

Co 

and  using  (17).  the  particular  solution  becomes 


My)  =  -f>  //- 

i=i  v 

1 


(x(T.x„.vn)|y) 


x/i(*n.Xn<  Va)dx„dv„.  (18) 


Assuming  now  an  initial  discretization  of  phase  space  (xj.  vk) 
and  applying  the  trapezoidal  rule,  the  integral  over  Ft  for  each 
species  /,  can  be  approximated  by 


N,  M, 

WjWk 

j=0  k=G 


G(x(r.Xj.vk)ly) 

fn 


j,  Vk) 


where  j  denotes  a  single  or  double  sum  approximating  a  single 
or  double  integral  in  space  and  k  denotes  a  single,  double,  or 
triple  sum  approximating  the  velocity  integral,  and  w/j,  are 
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Fig.  3.  Stretching  of  phase  space  and  panicle  insertion.  Black  dots:  Original 
panicles  White  dm:  Inserted  panicle. 

quadrature  weights.  Setting  7^  jk}  —  wj  Wk/i(^>*  xj *  Vk),  (18) 
is  approximated  by 


m  N,  M, 

<My}  =  -  EE  E  ,^i<J'klfi(x{T.xJ.v,)ly)  (19) 

J  j— 1>  k— 0 

where  7^  ,  j  k }  is  the  phase  space  weight  of  a  macro- particle*  The 
sum  in  (19)  resembles  the  sum  in  (14). 

Due  to  the  stretching  of  phase  space  during  a  simulation,  it 
often  happens  that  adjacent  Lagrangian  particles  on  a  curve  in 
phase  space  move  fur  apart,  as  depicted  in  Fig,  3.  In  this  case, 
ihe  quadrature  formula  becomes  inaccurate,  hut  this  can  be  con¬ 
trolled  by  inserting  new  particles  to  maintain  resolution.  Sup¬ 
pose  that  a  particle  needs  10  he  inserted  between  {i,  j,  k}  and 
{ /  ,  j  +  I ,  k),  based  on  some  criterion,  as  shown  in  Fig.  3.  Recall 
that  the  initial  phase  space  coordinates  are  Lagrangian  param¬ 
eters.  We  insert  a  new  particle  at  the  midpoint  of  the  interval 
between  {i,  j,  k}  and  {/',  j  +  I ,  k}  in  phase  space  at  t  =  L>.  The 
phase  space  weights  7  of  the  new  particle  and  the  existing  par¬ 
ticles  {7 ,  j ,  k(  and  {i,j  +  l.k}  are  determined  in  accord  with 
the  trapezoid  rule  over  the  initial  condition.  The  phase  space  co¬ 
ordinates  of  the  new  particle  at  time  T  are  determined  using  a 
cubic  interpolating  polynomial  with  respect  to  the  Lagrangian 
parameters  at  particles  { i,  j  —  L  k } ,  { 6  j .  k } .  { it  j  +  1 .  k  | .  and 
{  /,  j  +  2,  k}.  The  same  process  can  be  carried  out  to  refine  with 
respect  to  velocity  as  well.  It  should  be  noted  that*  in  principle, 
intermediate  times  could  be  used  to  define  the  Lagrangian  pa¬ 
rameter  for  point  insertion,  however,  lhat  is  not  done  here.  Fur¬ 
ther  details  about  the  particle  insertion  scheme  will  be  discussed 
elsewhere  [41  f  Related  work  is  in  references  [29]  and  [37], 

C.  Tree  code 

Evaluating  the  particular  solution  in  (14)  and  (19)  is  an  ex¬ 
ample  of  /V-body  problem  and  the  central  processing  unit  (CPU) 
time  is  an  important  issue  [42 1,  We  employ  a  treecode  to  reduce 
the  operation  count  from  0(N2)  to  ()(N  log  N).  In  a  treecode, 
the  particles  are  divided  into  a  hierarchy  of  clusters  and  the  par¬ 
ticle- particle  interactions  are  replaced  by  particle-cluster  inter¬ 
actions,  which  are  evaluated  using  multipole  expansions,  Barnes 
and  Hut  [261  used  monopole  approximations  and  a  divide-and- 
conquer  evaluation  strategy,  while  Greengard  and  Rokhiin  |27] 
used  higher  order  spherical  harmonics  expansions  and  a  more 
sophisticated  evaluation  procedure,  Trcecodes  have  been  very 
successful  in  particle  simulations  and  there  is  ongoing  interest 


Fig.  4,  Panicle-cluster  interaction:  panicle  y.  cluster  C.  cluster  center  x, 

in  optimizing  their  performance  [43]-[5l],  Our  approach  fol¬ 
lows  [29],  whieh  extends  the  Barnes-Hut  treecode  [26]  in  var¬ 
ious  ways. 

Part i cl e-Ch tste r  In terac ti ons:  The  part  ic u lar  sol u  1  i on  c an  be 
cast  in  the  form 

rn  N  t 

^/•(y)  —  -  E  — <My) =  E  wJG(xjiy)-  (2°) 

.=1  -n  3= 1 

The  simplest  procedure  for  evaluating  is  direct  summation, 

but  this  requires  ()(N2)  operations  which  is  prohibitively  ex¬ 
pensive  when  N  is  large.  The  key  idea  in  a  treecode  is  to  replace 
the  particle-particle  interactions  in  (20)  by  suitably  chosen  par- 
tic  le-elu$tcr  interactions. 

The  potential  <h,(y)  is  first  expressed  as 

$;(y)  =  EE"«Gt*iW  =  E*^-^  {2l) 

c  jec  c 

where 

$;{y>G)  =  E  wjG(xj|y)  (22) 

J€C 

is  the  potential  due  to  the  interaction  between  particle  y  and  a 
duster  of  panicles  C  —  {xjx,  £  C}  (see  Fig.  4)r  The  proce¬ 
dure  for  choosing  the  clusters  will  be  explained  below;  for  now, 
it  is  enough  to  assume  that  the  clusters  in  (21 )  are  nonoverlap¬ 
ping  and  their  union  is  the  entire  particle  distribution. 

Next,  perform  a  Taylor  expansion  of  the  Green's  function 
about  the  duster  center  xr 

$i{y,C)  =  E»j*(*  +  (*i  -  xc)|y) 

J€C 

~  E  E  TF^G(X-|y)(Xi  -  Xr)' 

j  €C  J= f) 

=  E  7fa*G(x'ly)  E  “  x'-)' 

/=0  ;6C 

P 

=  Y,Tl(xc,y)M,(C)  (23) 

f=0 

where  /ms  the  order  of  approximation,  7}  (x,. ,  y )  is  the  /  th  Taylor 
coefficient  of  the  Green's  function,  and  M\{C)  is  the  /th  mo¬ 
ment  of  the  duster.  Note  lhat  Cartesian  multi -index  notation 
is  used,  e.g.,  in  two  dimensions  if  we  set  x  —  (,T|.:r2)  and 
l  —  (/}  ^2),  then  xJ  =  and  so  forth.  The  error  in  the 

approximation  is  0((r/R)p )*  where  r  =  maxy  ||Xj  -  x,,||2  is 
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level  1 


level  2 


level  4 


Fig.  5,  Examples  of  iree  structure,  (a j  Standard  (h)  Adaptive  129]. 


function  ComputePotentialfy ,  C) 
is  the  Taylor  approximation  sufficiently  accurate? 
if  yes,  compute  $i{y,C)  by  Taylor  approximation 
if  no,  does  C  have  any  children? 

if  yes,  call  ComputePotentiaI{y.  O')  for  each 
child  C/  of  C 

if  no,  compute  $,(y,C)  hy  direct  summation 


Fig,  ft .  Proccd  u  rc  for  ev  a  I  uati  ng  a  pa  rti  ck-cl  us  ter  i  n  tcraction  •!* ,  ( y .  C ) 


the  duster  radius  and  R  —  J|xr  —  yj[s  is  the  particle-cluster 
distance.  The  speedup  occurs  because  the  cluster  moments  are 
independent  of  the  particle  y,  and  the  Taylor  coefficients  are  in¬ 
dependent  of  the  number  of  particles  in  the  cluster  C,  The  Taylor 
expansion  in  (23)  is  a  Cartesian  analog  of  the  classical  spherical 
harmonics  mullipole  expansion  [52].  A  recurrence  relation  can 
he  used  to  efficiently  compute  the  Taylor  coefficients  to  any  de¬ 
sired  order  [29]. 

Tree  Structure:  The  set  of  all  clusters  has  a  hierarchical  tree 
structure.  In  a  standard  treeeode,  the  clusters  on  each  level  of  the 
tree  are  uniform  cubes  obtained  from  bisecting  the  previous  gen¬ 
eration  of  clusters  in  each  coordinate  direction.  This  is  shown 
in  Fig.  5(a)  for  a  set  of  particles  on  a  spiral  curve.  However, 
in  many  cases,  it  is  beneficial  to  shrink  the  boxes  as  shown  in 
Fig.  5(h).  The  nonuni  form  clusters  in  the  second  scheme  are 
well -adapted  to  the  particle  distribution  and  this  can  lead  to  a 
savings  in  CPU  time. 

Evaluation  Strategy:  The  potential  is  evaluated  using 
a  recursive  divide-and -conquer  strategy  [26].  As  indicated  in 
(21),  the  potential  at  a  given  point  y  is  expressed  as  a  sum 
of  particle-cluster  interactions  C)  for  suitably  chosen 

clusters  C,  The  treecode  has  three  options  for  evaluating  each 
term  Taylor  approximation,  direct  summation,  or 

descending  the  tree.  The  present  simulations  used  the  standard 
multi  pole  acceptance  criterion  (MAC),  he.,  the  Taylor  expan¬ 
sion  is  used  only  if  r/R  <  0,  where  r  is  the  radius  of  the 
cluster,  R  is  the  distance  between  the  particle  and  the  cluster, 
and  0  is  a  user-specified  parameter  [26],  [45].  If  the  MAC  is  not 
satisfied,  the  code  considers  the  interactions  between  particle  y 
and  the  children  of  cluster  C.  If  the  cluster  C  has  no  children, 
then  it  is  a  leaf  of  the  tree  and  <E>z(y,  C)  is  evaluated  by  direct 
summation.  The  procedure  is  sketched  in  Fig.  6. 

Electrostatic  Force:  The  method  described  above  can  be 
adapted  to  compute  the  electrostatic  force  on  a  particle  x*  in 
time  0(N  log  iV),  Substituting  (23)  for  the  particle-cluster  in¬ 
teraction  into  (20)  for  the  particular  solution,  using  the  relation 

E.‘(x*J  =  -Vy<f>,(y)|y=Xt 

tor  the  electric  field  at  the  panicle  x*.,  and  noting  that 
VyO(x|y)  —  — VxG(x|y),  the  electric  field  is 

(VXfTt(xr.xk))Mt(C). 

C  l=o 

Then  the  electrostatic  force  on  x*  due  to  species  i  is  q^k E;, 
where  qk  Wk  is  the  total  charge  on  particle  k. 


Regularizing  the  Green's  Function:  In  l  wo  and  t  hree  di  men  - 
sions,  the  imerparticlc  force  F,  in  the  equations  of  motion  (2) 
becomes  singular  as  the  distance  between  the  particles  tends  to 
zero,  placing  a  constraint  on  the  maximum  allowable  time  step. 
Our  approach  to  overcoming  this  problem  is  to  regularize  the 
Green's  function  [29],  |37|,  In  two  dimensions,  we  use 

g.i»i»)--l°'(l|x7lli+<’) 

Hi T 

where  A  is  a  smoothing  parameter,  so  that  the  maximum  force 
is  proportional  to  1/A.  The  appropriate  value  of  A  could  be  de¬ 
termined  by  physical  properties  of  the  system,  e.g.,  in  a  neutral 
plasma,  one  could  use  the  Debye  shielding  length  to  choose  A. 
However,  in  the  work  presented  here,  A  was  treated  as  a  numer¬ 
ical  parameter 

/>.  Boundary  Integral/Homogeneous  Solution 

Having  completed  our  discussion  of  the  particular  solution, 
next  we  describe  a  boundary  integral  method  for  the  homoge¬ 
neous  solution;  for  an  alternative  finite  element  formulation,  see 
[38].  Recall  that  the  potential  function  4>  satisfies  the  Poisson 
equation  (4)  subject  to  Dirichlet,  Neumann,  or  mixed  boundary 
conditions  depending  on  the  application.  The  homogeneous  so¬ 
lution  is  given  in  ( 1 1 )  or  { 1 2)  and  at  each  point  on  the  boundary 
cither  V4>  ■  n,  or  *  n  is  specified.  As  an  ex¬ 

ample,  consider  the  case  of  Dirichlet  boundary  conditions  in 
which  —  nr  is  known  but  ■  n  =  0  is  unknown.  We  re¬ 
quire  the  boundary  condition  to  be  satisfied  in  the  sense  that 

Hm  4>(y)  =  n?(z),  for  z  €  Of];  (24) 

A  similar  relation  holds  as  well  for  the  other  possible  boundary 
conditions. 

For  our  Dirichlet  example,  this  yields  an  equation  for  0\  a  2 
x  2  linear  system  in  one  dimension 

lim  ^  0(zo)G(zo\y)) 

=  f-rrt(2i)  -  ^(zo))  -  lim  -  $r(y))  (25) 

\  *  l  )  V— *0,i> 

and  an  integral  equation  of  the  first  kind  in  higher  dimensions 

<f  (mC(x\z)<lSx  =  +  <&P{z) 

m 

+  /  rt(x)VxG(x|z)  ■  ih/.S’x  (2ft) 
on 

where  n;{z)  on  the  right  side  of  (26)  is  divided  by  2  to  account 
for  the  singular  nature  of  the  limit  as  y  approaches  the  boundary 
153]. 
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■  While  (25)  can  be  directly  solved  lor  the  unknowns.  (26)  will 
be  solved  by  a  collocation  method.  First,  discretize  the  boundary 
integrals  into  panels  for  j  =  1 . M 


M  ,  (  * 

£  /  P(*)GWi)d5»  =  -2SJ  +  *P( z) 


£  /  rt(x)VxG(x|z).nfiSx. 
=lsn 


(27) 


Choosing  z  “  z,  to  be  the  center  of  the  ah  panel  and  letting 
t\  —  a7  and  0  =  /Jj  be  constant  along  each  panel  #Qj  (for 
example),  we  obtain 


m  ,  , 

P<  Y1  ,"'g(xj'Iz.)  =  ~~~  +  M{z.) 

j=i  / 

M 

+  ^w,Vx(7(x,i|z,) -Iij,  (28) 

i=i  i 

where  Xj/  are  Gaussian  quadrature  points  on  the  jth  panel  and 
wi  are  the  corresponding  weights.  This  yields  a  linear  system 
for  the  0] 

A0  =  k 

where 

<’ij  =  ^w(6'(x,i|z,) 

; 

and 


i  h.  J XI 

K|  =  —~2  +  *Mzi)  +  ^«j^-*«|VxG{xj(.|z<)  ■  n_j. 

j=i  i 

This  formulation  worked  well  in  our  preliminary  work  [36]. 
For  geometrically  complex  domains,  the  linear  system  becomes 
challenging  to  solve  and  an  efficient  method  such  as  precondi¬ 
tioned  generalized  minimal  residual  (GMRES)  is  required.  Due 
to  the  special  form  of  the  matrix  elements  atJ,  the  matrix -vector 
multiplication  in  each  step  of  GMRES  can  be  computed  effi¬ 
ciently  using  the  treecode.  Hence,  the  matrix  A  never  needs  to 
he  explicitly  constructed  and  stored,  making  geometrically  com¬ 
plex  domains  more  tractable. 

This  approach  can  also  be  used  for  problems  with  Dirichlet 
and  Neumann  conditions  on  different  portions  of  the  boundary 
f+)S7  =  dUo  U  flfijv.  In  this  case,  we  represent  <!>//  as  a  combi¬ 
nation  of  single  and  double  layer  potentials 

<My)=  f  rv(x)G{x|y)rfSx  +  f  (i(x)YxG(x\y)-ndSx 

DR.v  onR 

where  f>,  /?  are  unknown  functions.  This  eliminates  the  need 
to  compute  a  boundary  integral  for  each  We  found  through 
numerical  experiments  that  the  resulting  matrix  A  is  well-con¬ 
ditioned  when  a  double  layer  potential  is  used  for  the  Dirichlet 


boundary  and  a  single  layer  potential  is  used  for  the  Neumann 
boundary.  Mixed  boundary  conditions  can  be  handled  by  going 
hack  to  the  original  formulation  using  ( 1 2)  for  $//. 


in.  Applications 

In  this  section,  we  give  an  overview  of  the  problems  to  which 
the  BIT  approach  has  been  applied.  We  divide  the  discussion 
into  two  sections  based  on  the  problem  dimension. 


A,  / -D  Problems 

We  used  BIT  to  simulate  several  benchmark  1-D  problems 
in  plasma  physics  including  sheath  formation,  virtual  cathode 
formation,  and  ihe  two- stream  instability.  In  the  latter  case,  we 
also  used  particle  insertion. 

Sheath  Formation:  The  1-D  sheath  domain  is  bounded  by 
two  parallel  metal  plates  a  distance  Z\  -  zQ  “  L  apart,  each 
held  at  a  potential  of  zero  volts.  The  domain  initially  contains  a 
neutral  plasma  with  equal  electron  and  ion  densities.  The  elec¬ 
trons  are  lighter  and  more  mobile  than  the  ions,  and  a  fraction  of 
the  electrons  near  the  wall  impact  the  wall  on  a  timescale  much 
faster  than  the  ion  transit  lime.  The  reduced  electron  density  near 
the  wall  leaves  a  net  positive  charge  in  a  n arrow  region  between 
the  wall  and  the  main  body  of  the  plasma.  This  gives  rise  to  a  po¬ 
tential  barrier  (sheath)  at  the  walls,  and  if  any  of  the  remaining 
electrons  are  to  escape,  they  must  have  a  large  enough  velocity 
to  overcome  the  sheath  potential.  Eventually  the  system  reaches 
a  steady  state  in  which  the  electrons  are  electrostatically  con¬ 
fined  to  the  interior  of  the  domain. 

Wc  looked  at  this  application  in  order  to  demonstrate  that  BIT 
can  reporoduce  the  undriven  sheath  in  a  decaying  plasma  as  an 
initial  benchmark.  This  sheath  decay  includes  some  kinetic  ele¬ 
ments,  since  the  tail  is  lost  first  and  the  Coulomb  collision  rate  at 
these  lowr  densities  is  insufficient  to  maintain  a  Maxwdl-Boll/- 
mann  distribution,  so  kinetic  effects  are  important,  and  hence 
this  is  best  represented  by  akinetic  model.  This  type  of  problem, 
decay  of  an  undriven  bounded  plasma  and  resulting  sheath  for¬ 
mation,  has  been  studied  by  PTC  at  some  length  [54H56J,  and 
is  widely  used  as  a  benchmark  for  other  models  including  fluid 
models  [57],  [58]. 

In  this  simulation,  the  electrons  are  represented  as  mobile  par¬ 
ticles,  while  the  ions,  being  much  heavier,  are  immobile.  The 
Poisson  equation  lakes  the  form 


vH  =  - 


N. 

-  x})  -  «, 
j=i 


(29) 


where  r.  is  the  electron  charge,  Nr  is  the  numher  of  electrons 
in  the  simulation,  and  nx  is  the  constant  background  ion  den¬ 
sity,  chosen  such  that  niL  —  w}Nr.  The  electrons  are  initially 
distributed  uniformly  in  space  with  a  random  thermal  velocity 
sampled  from  a  Maxwellian  distribution. 

The  potential  is  the  sum  of  three  terms 


${?/}  =  ‘Mi/)  +  ‘My}  +  M(y) 


(30) 
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Position  (m) 


Fig,  7.  Potential  t  as  a  function  of  position  across  domain  for  BIT  and  PIC. 


Fig,  8,  Electron  density  for  BIT  and  PIC,  Ion  density  is  also  shown. 


where  is  due  to  the  electrons,  is  due  to  the  background  ion 
density,  and  is  the  homogeneous  solution  needed  to  enforce 
the  boundary  conditions 


o 

-[*«  +  *;]{*(,.*, )  - 


(31) 


The  term  is  given  by  (14}T  <1 >,  is  obtained  by  integrating 
twice,  and  $//  is  found  by  substituting  (25)  into  (t  1) 


*.(»)  =  ) 

f  (32i 

•*(»)  =  C(*l»W(*o)-C(*ilirW(*i)  ) 


where  we  have  also  imposed  the  boundary  conditions  <b(20)  — 
fb(*i)  —  0  in  (J  l).  With  these  boundary  conditions,  the  limits 
in  (25)  give 


/#(*„)  =  -Gfzolz,)-1  [M*i)  +  *,(*,)]  1 
0{zi)  ~  G(zil^a)  1  ['M^n)  +  ^t(^o)]  J 

BIT  results  will  be  compared  with  PIC  results  from  [59],  Note 
that  the  l-D  Green's  function  in  (7)  is  piecewise  linear  and  so 
the  first  order  Taylor  approximation  is  exact;  this  implies  that 
BIT  computes  exact  forces  and  potentials  in  one  dimension, 
and  gives  identical  results  to  direct  summation.  In  the  work  pre¬ 
sented  here,  the  BIT  simulation  was  initialized  using  the  same 
parameter  set  as  in  the  PIC  simulation  so  that  the  only  difference 
Wits  in  the  field  solver  In  both  simulations,  the  domain  length 
is  L  -  0.1  in,  the  number  density  of  each  species  is  1013  m”3, 
and  the  electron  temperature  is  l  eV.  The  initial  number  of  elec¬ 
trons  is  4000  in  both  simulations. 

Fig.  7  plots  the  potential  as  a  function  of  position  across 
the  domain  for  BIT  and  PIC.  There  is  no  significant  difference 
between  the  two  methods,  either  in  magnitude  or  sheath  thick¬ 
ness.  The  calculated  matrix  sheath  thickness  for  this  domain  is 
5.4  mm  [60],  comparable  to  the  computed  results  of  approxi¬ 
mately  10  mm  for  both  BIT  and  PIC, 


Fig*  9*  Number  of  electrons  as  a  fund  ion  of  time  for  BIT  and  PIC 


Fig.  8  plots  the  electron  and  ion  densities  for  BIT  and  PIC. 
As  for  the  potential,  the  computed  particle  densities  for  the  two 
methods  are  almost  the  same* 

Fig.  9  plots  the  number  of  electrons  in  the  simulation  as  a 
function  of  time.  Again  the  agreement  between  the  BIT  and 
PIC  results  is  very  good.  Both  have  an  initial  depletion  rate 
that  matches  the  theoretical  result.  The  particle  number  drops 
slightly  faster  for  BIT,  resulting  in  approximately  20  fewer 
particles  after  convergence*  hut  this  has  negligible  effect  on 
the  potential  and  the  panicle  densities*  To  examine  whether 
the  additional  loss  in  BIT  arises  from  the  initial  seed  of  the 
Maxwellian  velocity  distribution,  many  additional  runs  were 
done*  each  using  a  different  seed  in  the  random  number  gen¬ 
erator.  All  of  the  results  were  within  1-3  particles  of  the  BIT 
results  in  Fig.  9,  Note  that  this  does  not  imply  that  BIT  and  PIC 
will  agree  in  a  full  discharge  simulation.  One  explanation  for 
the  discrepancy  in  the  number  of  particles  is  Lhai  BIT  includes 
short  range  Coulomb  collisions  which  are  neglected  in  PIC* 
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Hence,  we  expect  that  ihe  tail  in  BIT  will  he  rcpupulaled  mure 
quickly  than  in  PIC.  Since  the  tail  is  the  only  pan  that  can 
escape  the  potential  well,  wc  then  expect  BIT  to  have  a  higher 
loss  rate  and  hence  a  slightly  lower  number  of  particles  than 
PIC. 

Virtual  Cathode:  Virtual  cathodes  (VC)  are  basic  structures 
that  arise  in  many  plasma  applications,  e.g.,  they  have  been  used 
to  explain  anomalous  acceleration  of  positive  ions  in  a  vacuum 
diode  [61],  and  have  also  been  proposed  as  a  confinement  mech¬ 
anism  for  controlled  fusion  [63].  The  formation  of  a  VC  is  un¬ 
derstood  in  general  terms,  but  predicting  the  details  is  difficult 
because  the  process  is  highly  nonlinear.  Such  systems  are  often 
simulated  using  PIC  (64],  but  it  is  not  clear  that  this  is  optimal 
since  sharp  gradients  in  the  plasma  properties  occur  near  the 
VC. 

In  previous  work,  wc  simulated  the  formation  of  a  virtual 
cathode  in  one  dimension  using  a  BIT  algorithm  [35).  Two  infi¬ 
nite  parallel  flat  plates  have  applied  voltages.  The  gap  between 
the  plates  is  initially  a  vacuum,  and  one  of  the  plates  is  healed  to 
Ihe  point  where  electrons  are  emitted  and  a  current  flows  across 
the  gap.  For  a  sufficiently  high  emission  rate,  the  maximum  cur¬ 
rent  density  between  the  plates  approaches  a  limiting  value  -/cl 
predicted  by  the  Child-Langmuir  taw.  When  this  happens,  a  vir¬ 
tual  cathode  is  formed  somewhere  between  the  two  plates,  and 
the  emitted  electrons  turn  around  and  strike  the  emitting  plate. 
As  the  electrons  repel  one  another,  the  potential  minimum  rises, 
thereby  allowing  the  emitted  electrons  to  again  start  crossing  the 
gap,  and  the  entire  sequence  repeats.  Wc  considered  the  classical 
formulation  in  which  the  electrodes  are  held  at  ground  and  the 
emitted  electrons  have  an  initial  velocity  [39). 

For  this  nonlinear  oscillator,  we  compared  PIC,  BIT,  and  di¬ 
rect  summation  (DS).  As  noted  in  the  section  on  sheath  forma¬ 
tion,  DS  and  BIT  agree  to  machine  round  off,  so  the  main  reason 
for  the  DS  calculation  is  for  timing  comparisons.  For  the  BIT 
and  DS  runs,  the  grid-free  solution  to  the  Poisson's  equation  is 
identical  to  (30),  except  that  nt  is  identically  zero.  The  PIC  rou¬ 
tine  employed  here  used  standard  area-weighting.  The  injected 
current  was  varied  by  adjusting  the  charge/mass  ratio  of  the  par¬ 
ticles.  With  a  small  injected  current,  the  solution  converged  to  a 
steady  state  and  all  three  methods  gave  the  same  result  to  within 
plotting  accuracy.  The  mesh  size  for  the  PIC  simulation  was  de¬ 
termined  by  decreasing  the  value  until  the  solution  converged 
to  the  nonoscillaiory  steady  state.  This  mesh  spacing  was  then 
used  for  the  remainder  of  the  PIC  simulations. 

With  a  higher  injected  current,  a  VC  forms  and  the  solution 
converges  Lo  a  time-periodic  state  with  nonuni  lorm  spalial  struc¬ 
ture.  Fig.  10  show*  the  time  trace  of  the  particle  density  {nfno) 
at  three  locations  xjL  between  the  plates;  just  past  injection, 
near  the  VC,  and  downstream  of  the  VC.  Fig.  10  reveals  signifi¬ 
cant  differences  between  the  BIT/DS  and  PIC  results.  First,  the 
oscillation  period  is  59c  greater  for  PIC  than  for  BIT/DS,  pro¬ 
ducing  a  significant  phase  lag.  Fig.  11,  shows  the  particle  distri¬ 
bution  in  phase  space  at  a  given  time.  The  plot  shows  that  some 
of  the  emitted  particles  turn  hack  and  approach  the  emitter  plate. 
The  PIC  results  have  irregular  voids,  while  the  BIT/DS  results 
have  a  smooth  particle  distribution. 

All  three  methods  used  the  same  particle  push  and  particle 
injection/deletion  procedures.  The  only  difference  between  the 
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time  t/T 

Fig.  10.  Time  trace  in  the  formation  of  a  virtual  cathode  ai  three  spatial 
locations  across  the  channel,  r/L  —  1/16. 5/1G.9/1G  Solid  line  is  RH7DS 
and  the  dashed  line  is  PIC. 


BIT 


Fig.  1 1 ,  Phase  space  plot  at  time  t/T  —  U.5.  Irregular  voids  appear  for  PIC. 
while  HIT/DS  has  a  smooth  panicle  distribution. 


three  routines  was  the  force  evaluation  kernel.  In  all  simulation 
results  presented,  the  time  step  used  was  that  of  the  time  con¬ 
verged  DS  method.  The  criterion  for  the  lime  step  was  that  the 
phase  space  point  wise  L2  norm  must  be  less  then  a  user  set  tol¬ 
erance 


DS{At)-Ds(^j 


<  0 


at  any  observed  time  between  tn  and  //.  This  guarantees  point- 
wise  convergence  in  phase  space.  The  observed  differences  in 
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Fig.  1 2.  Two-stream  instability:  initial  electron  distribution  in  phase  space. 


Lht;  simulation  results,  shown  in  Figs.  10  and  11,  can  be  at¬ 
tributed  Lo  the  error  in  intcrparticlc  force  within  a  mesh  cell  of 
a  PIC  simulation,  as  shown  in  Fig,  1.  These  preliminary  runs 
were  performed  using  Matlab;  PIC  required  7  h,  BIT  12  h,  and 
DS  seven  days.  Further  work  is  needed  to  optimize  the  run  time 
of  the  BIT  code,  but  the  results  suggest  that  BIT  is  capable  of  re¬ 
solving  small-scale  features  better  than  a  mesh-based  PIC  code, 

Two-Stream  Instability:  This  example  concerns  a  neutral 
plasma  with  equal  numbers  of  ions  and  electrons  in  a  1-D 
periodic  domain.  As  above,  the  electrons  evolve  dynamically 
and  the  ions  are  immobile.  Consider  an  electron  distribution 
function  of  the  form  f(x,  w,fo)  =  f(x.  — t;,  to)  =  j/flw]}.  i.e., 
independent  of  position  x  and  an  even  function  of  velocity  ik 
This  is  an  unstable  equilibrium  solution  of  the  Vlasov-Poisson 
system,  A  perturbation  in  velocity  gives  the  electrons  a  spatial 
variation  in  momentum.  This  leads  to  nonuni formity  in  the 
electron  spatial  distribution,  thereby  generating  an  electric  field 
which  feeds  back  into  the  velocity  perturbation,  causing  phase 
space  to  fold  up  on  itself. 

As  with  the  virtual  cathode,  we  focus  on  a  cold  two-stream 
instability  in  the  sense  thai 


s(M)  = 


if  v  = 
otherwise 


where  v„  is  the  initial  unperturbed  speed  of  the  electrons.  The 
initial  velocity  perturbation  has  the  form 

f  \  •  (*% 

—  asm  f  Ztt - 

V  -l  —  2q 


where  a  is  the  perturbation  amplitude.  Fig,  12  shows  the  initial 
electron  distribution  in  phase  space.  The  electrons  are  modeled 
as  panicles  while  the  ions  are  treated  as  a  constant  background 
density  ti.,\  The  potential  of  the  domain  is  then  described  by 


V2 


j=i 


with  periodic  boundary  conditions 


dx 


dx 


where  r  is  the  charge  on  an  electron  and  Nr  is  the  number  of 
electrons  in  the  simulation. 

As  in  the  sheath  formation  problem,  the  particular  solution 
4>p  is  the  sum  of  <f>r  and  from  (32),  The  boundary  condi¬ 
tions  determine  the  homogenous  solution  given  by  (25),  Le>, 
*(*)  =  ${21)*  d$(z0)/dT  =  d${zt)fd,r  gives  = 

n{zx),  ft(zo)  -  (}(zi).  Thus,  (25)  gives 


0(*i)C(*i|*o)  =  -  |«(*))  + 

and 

Using  the  periodic  boundary  conditions  and  noting  that 
=  G'(20|*i).  we  obtain 

™(*i)  =  “  ($p(zo)  +  (hr(^i )) 


and 


PM 


i 

2G(z\  |^o ) 


(*p(*o)  ~  $p(*i))  ■ 


The  initial  ion  and  electron  density  is  4  x  10s  cnr3.  All  the 
phase  space  plots  are  normalized  to  the  initial  electron  velocity 
?/  =  c.  Time  stepping  was  performed  using  the  fourth -order 
Runge-Kuita  (RK4)  method  and  the  results  were  considered 
time  accurate  when  the  difference  between  runs  with  A/  and 
Af/2  was  less  then  5%.  Spatial  (i.e.,  phase  space)  convergence 
was  achieved  by  systematically  increasing  the  number  of  sim¬ 
ulation  points  Nr.  Fig.  13,  lop  frame,  shows  a  time  converged 
result  for  Nr  =  3200  after  64000  time  steps  with  A t  =  (r  — 
zo)/kc<  where  k  —  Hi  x  10,  The  electron  distribution  function 
is  concentrated  on  a  rollcd-up  curve  in  phase  space.  Fig.  13, 
middle  and  bottom  frame,  shows  that  the  simulation  loses  spa¬ 
tial  accuracy  later  in  time  when  adjacent  particles  become  sep¬ 
arated. 

To  overcome  the  loss  of  resolution  that  inevitably  occurs 
using  a  fixed  number  of  particles,  we  implemented  a  particle 
insertion  scheme  as  discussed  briefly  in  Section  II-B.  In  these 
results  a  simple  distance  criterion  is  used,  i.e.,  a  new  particle  is 
inserted  between  particles  i  and  i  4  I  when 


y/(Xi  -  Xi  +  i)2  4  (in  -  vl  +  l)2  >  -  S 

where  S  is  the  initial  panicle  spacing.  Fig.  14  shows  the  re¬ 
sult  obtained  using  particle  insertion,  starting  wiih  200  points 
and  using  the  same  time  step  as  in  Fig,  1 3,  This  simulation  was 
checked  for  convergence  both  in  time  as  well  as  in  space  by 
computing  with  twice  as  many  iniiial  points.  The  results  indi¬ 
cate  that  particle  insertion  is  an  efficient  means  for  maintaining 
accuracy  as  the  solution  evolves  inti  me.  A  detailed  analysis  of 
these  results,  as  w'dl  as  extension  to  the  warm  two-stream  insta¬ 
bility,  arc  topics  for  future  work. 
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3200  points  t=2.5s 
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3200  points  t=3.3s 
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Fi  g.  I ,T  TWo-  sire  am  i  nsiabili  ly :  phase  s  pa  cc  plot  wi  th  ;Vr  =  3201) , Top  frame 
t  —  2,5,  Middle  frame  *  =  3  3,  Bottom  frame  f  —  4.3, 


B,  2-D  Problems 

In  addition  to  the  l-D  problems  discussed  above,  we  also 
applied  the  BIT  approach  to  several  2-D  problems,  including 
planar  and  cylindrical  ion  optics,  and  particle  dynamics  in  a  Pen- 
ning-Malmberg  trap.  We  performed  a  detailed  timing  compar¬ 
ison  in  2-D  planar  geometry  between  BIT  and  multigrid  PIC, 
We  start  with  the  timing  comparison  and  then  discuss  the  appli¬ 
cations  listed  above, 

77m mg  Resi dts;  In  this  s  ection,  we  c  om  pare  thee  f fi  c  ie  nc  y  o  f 
BIT  and  PIC  simulations  for  a  2-D  test  case,  focusing  on  the 
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Fig,  14,  Two-stream  instability;  phase  space  plot  using  point  insertion,  initial 
Nr  —  200.  Top  frame  t  =  2.5.  Butin  m  frame  t  =  3,3. 


question:  for  a  given  level  of  accuracy,  which  method  is  faster? 
In  this  test,  N  particles  were  randomly  placed  in  a  square  do¬ 
main  with  grounded  boundary  and  the  electric  field  was  com¬ 
puted  at  each  panicle  location  using  BIT,  PIC,  and  direct  sum. 

Note  that  if  the  homogeneous  solution  in  (12)  is  computed 
exactly,  then  direct  sum  yields  the  exact  solution  of  the  test 
problem.  The  Green's  function  for  a  square  domain  can  be 
expressed  as  an  infinite  series,  but  we  chose  to  satisfy  the 
boundary  condition  using  the  frce-space  Green's  function.  This 
requires  solving  a  boundary  integral  equation  for  the  homoge¬ 
neous  solution.  We  expressed  the  homogeneous  solution  as  a 
single  layer  potential  and  subdivided  the  boundary  into  panels 
as  discussed  in  Section  II-D.  In  this  example,  the  terms  in  the 
matrix,  at}  -  [  G{xi\x3)dSrr  can  he  integrated  exactly  and 
the  panel  strengths  are  determined  to  within  roundoff  error.  The 
number  of  panels  in  the  direct  sum  solution  is  chosen  so  that 
the  error  is  less  than  1  %.  The  resulting  boundary  integral  direct 
sum  result  is  taken  as  die  exact  solution. 

The  error  is  measured  by  the  expression 


error 


t  jY 


E'^  -  £" 


1=1 


where  Ef*  is  the  direct  sum  electric  held  at  particle  i  and  E-wr 
is  the  BIT  or  PIC  electric  field  approximation.  Table  I  shows 
liming  results  versus  number  of  panicles  for  D55,  PIC- MG, 
and  BIT.  PIC-M  denotes  a  PIC  solution  on  an  M  x  M  grid 
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TABLE  [ 

CPU  Time  (a)  i-tm  DS.  PIC-M.  B1T-T  (See  Text  11m  Dhi  initiuns) 


N 

DS 

P1C-12B 

P1C-2S6 

PIC-512 

HIT-1 

BIT-2 

150K 

792 

0.50 

4.53 

32.4 

19 

22 

JOCK 

2995 

0,52 

4.54 

32.7 

38 

45 

6{X1K 

ii  m 

0.55 

4 .55 

32.9 

78 

90 

TABLE  11 

Error  tor  PIC-M,  BIT-T  (Siu;.  TEXT  lor  Ditinhions) 


N 

DS 

PEC-1 2K 

PIC-256 

PIC-512 

BIT- 1 

BIT-2 

150K 

- 

m 

10% 

77c 

8% 

“<  1% 

300K 

- 

13% 

m 

6% 

7% 

<  1% 

600K 

- 

14% 

10% 

8% 

6% 

<  1% 

using  a  multigrid  solver  [62].  BIT-T  denotes  a  boundary 
integral/ircecode  solution  using  a  Taylor  expansion  of  order 
T  (BIT-1  used  MAC  parameter  ft  —  0,9  and  BIT-2  used 
#  -  0.45).  Table  II  shows  the  corresponding  errors. 

A  few  observations  can  be  made.  1)  Asa  function  of  N ,  DS 
is  ( ){N 2),  PIC-MG  is  almost  independeni  of  N  over  this  range, 
and  BIT  displays  O(N)  behavior  although  theoretically  it  is 
()(N  log  N).  2)  For  fixed  jV,  the  time  required  for  PIC  increases 
by  a  factor  of  8  when  the  mesh  is  refined  by  a  factor  of  2  in  each 
direction.  Alternatively,  the  lime  required  fur  BIT  increases  by 
a  fraction  less  than  one  when  going  from  first  to  second  order 
Taylor  expansion.  In  terms  of  accuracy*  ihe  PIC-512  and  BIT- 1 
results  have  comparable  errors,  while  the  BIT-2  results  have 
much  smaller  errors.  Note  that  PIC  would  require  an  exceed¬ 
ingly  line  mesh  to  achieve  the  same  accuracy  as  BIT-2,  with  a 
correspondingly  large  increase  in  the  CPU  time  (estimated  lo 
be  ai  least  2000  s).  Note  that  in  Table  II,  ihe  PIC  results  display 
firsi-order  convergence;  this  is  due  lo  the  choice  of  interpolation 
and  ihe  effect  of  differencing  the  potential  to  obtain  the  electric 
field. 

The  timing  results  imply  that  for  problems  where  strong 
coulomb  interactions  matter,  an  extremely  fine  spatial  dis¬ 
cretization  may  be  required  to  resolve  these  effects  in  a  PIC 
simulation.  In  summary,  for  the  choice  of  parameters  used  in 
this  strong  coulomb  interaction  test,  PIC  is  more  efficient  for 
10%  accuracy  and  BIT  is  more  efficient  for  1%  accuracy.  This 
implies  that  for  a  large  class  of  problems,  where  sharp  gradients 
in  plasma  properties  develop  in  localized  regions,  it  may  be 
ideal  to  think  about  a  hybrid  PIC-BIT  algorithm.  In  such  an 
algorithm,  BIT  would  only  be  applied  to  small  patches  of  the 
simulation  domain  in  order  to  correct  the  PIC  calculation. 

Planar  Optics:  Ion  thrusters  for  spacecraft  propulsion  op¬ 
erate  by  electrostatically  accelerating  ions  through  a  set  of  ion 
optics.  The  optics  consist  of  a  screen  grid  and  an  accelerator 
grid,  each  of  which  typically  has  many  thousands  of  hexago¬ 
nal  ly-arranged  apertures.  Ion  optics  simulations  usually  focus 
on  a  single  aperture  and  the  immediate  vicinity  upstream  and 
downstream.  Fig.  15  show’s  the  geometry  and  boundary  condi¬ 
tions  in  the  present  example.  Ions  are  introduced  at  Ihe  upstream 
edge  of  the  domain  at  a  discharge  potential  of  about  2000  V.  The 
screen  grid  has  a  potential  25  V  below  the  discharge  potential 
and  the  accelerator  grid  has  a  negative  potential  on  the  order 
of  -200  V.  Ions  are  accelerated  electrostatically  by  the  drop  in 
potential,  providing  thrust.  Downstream  of  the  optics,  the  beam 
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Fig  15  Geometry  for  planar  ion  optics  simulation.  D'  Pmchlet.  N  Neumann, 
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Fig,  16  Snapshot  of  particles  in  a  planar  ion  optics  simulation  Also  shown 
are  ihe  leaf  clusiers  in  the  treecode,  assuming  a  maximum  of  2(XI  particles  per 
leaf. 


is  neutralized  by  an  electron  source,  making  the  downstream 
potential  approximately  zero,  although  ihe  present  BIT  simula¬ 
tions  included  only  the  ions. 

In  the  present  simulations,  ions  were  injected  at  a  rale  of 
10  per  iteration,  each  with  a  weight  corresponding  to  number 
density  2.3  -  10 13  in'3.  This  low  density  was  used  because 
the  absence  of  neutralizing  electrons  causes  the  space -charge 
limit  to  he  reached  at  much  lower  density  than  in  typical  ion 
thruster  operation.  The  timestep  w'as  At  =  1()-9  s  and  the  par¬ 
ticles  were  injected  with  Maxwellian  velocities  centered  about 
2000  m7s.  The  boundaries  were  modeled  using  a  combination 
of  single  and  double  layer  potentials,  with  panels  uniformly  dis¬ 
tributed  on  each  segment.  The  line  integrals  of  the  Green's  func¬ 
tion  and  its  normal  derivative  were  integrated  exactly  over  each 
panel — our  previous  work  used  quadrature  [36],  but  exact  inte¬ 
gration  was  found  to  mitigate  the  effect  of  the  kernel  singular¬ 
ities  in  evaluating  particle- boundary  interactions.  The  treecode 
used  a  fourth -order  Taylor  expansion  and  the  MAC  parameter 
was  0  =  0.1. 

Fig.  16  shows  a  snapshot  of  the  computed  particle  locations. 
Also  shown  arc  the  leaf  clusters  in  the  treecode,  assuming  a 
maximum  of  200  particles  per  leaf.  The  clusters  conform  to  the 
particle  locations:  no  effort  is  expended  on  empty  regions,  re¬ 
gions  with  low  particle  density  have  fewer  clusters,  and  regions 
with  high  particle  density  have  many  clusters,  A  detailed  anal¬ 
ysis  of  the  results  is  in  preparation. 

Cylindrical  Optics:  In  this  case,  BIT  is  applied  to  a  mure 
complex  problem  than  in  the  previous  test  cases.  One  problem 
that  arises  is  the  inclusion  of  an  electron  population  for  the  direct 
summation  and  BIT.  The  PIC  potential  solver  models  a  Boltz¬ 
mann  electron  fluid  at  the  mesh  nodes,  but  BIT  has  no  mesh  for 
such  a  simulation  and  there  is  no  known  Green's  function  for 
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the  fully  nonlinear  problem.  The  free- space  Green's  function 
for  the  axisym metric  case  is 


c(x'y)=i/ 


dB 

(i  -  m  sin2  0)4 


Kjm) 
7 \\fL 


(34) 


where  x,  y  have  cylindrical  coordinates  (rx, 2X),  (ru,  zu),  L2  - 
(rif  +  rv)2H-(2x-Z|f)a;m  =  4r,rtf/La»and/f(m)  is  the  com¬ 
plete  elliptic  integral  of  the  first  kind.  The  present  axi symmetric 
treccodc  is  relatively  inefficient  compared  to  the  planar  version 
due  to  the  tack  of  a  recurrence  relation  for  the  derivatives  of 
the  axisymmetric  Green's  function.  As  a  result,  in  this  example, 
the  treecode  CPU  lime  is  at  least  an  order  of  magnitude  greater 
than  the  PIC  CPU  time.  The  treccodc  could  have  been  made 
more  efficient  by  using  a  lookup  table  for  the  elliptic  function 
evaluations  or  a  different  form  of  the  Green's  function  |65|,  The 
problem  does  not  arise  in  three  dimensions  because  an  efficient 
recurrence  relation  is  known  for  that  case  [29] * 

The  simulation  domain  radius  was  10” 3  in,  the  upstream 
length  was  2  x  10“3  mt  the  downstream  length  was  4  x  10”3  m, 
Lhe  screen  grid  was  0,4  mm  thick  with  a  i  .6-mm-diameler  aper¬ 
ture,  and  the  aecel  grid  was  0.8  mm  thick  with  a  1.0-mm-diam- 
clcr  aperture.  As  for  planar  optics,  the  boundary  was  modeled 
as  a  combination  of  single  and  double  layer  potentials.  The 
boundary  matrix  A  was  constructed  using  an  eight  point 
Gaussian  quadrature  rule.  To  control  errors  at  sharp  comers,  a 
Chebyshev  panel  spacing  was  used  along  boundary  segments. 
The  total  number  of  panels  was  512.  The  upstream  domain 
boundary  potential  was  1800  V,  the  screen  grid  was  1775  V, 
the  acccl  grid  was  -210  V,  the  downstream  boundary  was 
22  V  (all  Dirichlet),  and  the  other  upper  boundaries  were 
homogeneous  Neumann,  The  smoothing  parameter  value  was 
S  =  4  x  HU5  m*  The  treecode  used  fourth  order  Taylor 
expansions,  the  MAC  parameter  was  B  —  CLIO,  and  the  leaves 
had  no  more  than  eight  particles  per  lowest  level  duster. 

The  comparison  done  here  is  between  direct  summation  and 
PIC  in  order  to  provide  a  baseline  comparison*  Because  the  dif¬ 
ferences  are  fairly  large  in  this  case,  a  comparison  using  the 
treccodc  would  appear  exactly  the  same  due  to  its  closeness  to 
direct  summation.  The  direct  summation  comparison  to  PIC  im¬ 
parts  all  the  relevant  information* 

Fig.  17  plots  the  average  relative  difference  in  the  particle 
force  between  PIC  and  direct  summation  and  the  average  elec¬ 
tric  field  across  the  domain,  In  the  central  region,  the  imposed 
electric  field  dominates,  giving  a  small  relative  difference.  In  the 
neutral  regions  of  the  domain  upstream  and  downstream  of  the 
ion  optics,  the  electric  field  is  small  on  average.  PIC  computes 
very  small  forces  in  these  regions,  while  direct  summation  still 
sees  large  interparticle  forces*  The  average  PIC  force  magnitude 
in  the  upstream  region  is  on  the  order  of  10“ 1 7  N,  while  the  av¬ 
erage  direct  sum  force  is  approximately  10“ 16  N.  A  detailed 
analysis  of  the  results  is  in  preparation. 

Penning-Matmberg  Trap:  A  Penning-Malmberg  trap  is  a 
grounded  conducting  cylinder  used  for  confining  an  electron 
plasma  1 66],  [67 1.  The  electrons  are  confined  by  applying  a 
magnetic  field  along  the  cylinder  axis  and  holding  the  cylinder 
end  caps  at  constant  voltage,  so  that  the  electrons  bounce 
back  and  forth  along  the  magnetic  field  lines  {see  Fig.  18). 
For  example  in  [67],  the  end  caps  were  held  at  50  V  and  the 
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Fig.  17.  Relative  z  force  difference  between  PJC  and  direct  summation, 


Fig  18  Sketch  of  a  Penning-Malmberg  trap  for  confining  an  electron  plasma, 

magnetic  field  strength  was  B  1  T*  Under  these  conditions, 
the  time  required  for  an  electron  to  complete  one  bounce  is 
much  smaller  than  the  characteristic  E  x  B  time  scale.  In  this 
case,  the  plasma  is  well  described  by  a  2-D  particle  model  in 
which  the  electrons  behave  like  line  charges  being  conveeted 
with  velocity  v  =  E  x  B/B 2  [68|.  The  system  has  been 
investigated  experimentally  and  computationally,  and  many 
interesting  phenomena  were  revealed  including  metastable 
crystalline  states  and  complex  dynamics  [66],  [67]* 

In  our  simulations,  the  applied  magnetic  field  B  is  a  specified 
constant  and  the  electric  field  E  is  computed  using  BIT.  In  addi¬ 
tion,  the  particle  insertion  scheme  discussed  in  Section  II- B  was 
used  to  maintain  resolution  of  the  electron  density.  The  insertion 
criterion  compares  the  distance  between  hypothetical  particles 
inserted  using  linear  interpolation  and  cubic  interpolation.  If  the 
distance  is  greater  than  5%  of  the  initial  particle  spacing,  a  new 
particle  is  inserted.  The  time  integration  method  used  is  RK4, 
The  BIT  field  evaluation  used  an  eighth  order  Taylor  expansion 
and  the  MAC  parameter  was  B  =  0.5.  The  Dirich  let  boundary 
condition  on  the  cylinder  wall  was  imposed  using  a  single  layer 
potential  and  the  terms  in  the  boundary  matrix  A  were  evaluated 
using  eighi-point  Gaussian  quadrature.  Refinement  tests  were 
carried  out  in  space  and  time. 

Here,  we  present  a  simulation  of  wave  breaking  in  which  the 
cylinder  parameters,  .size,  magnetic  field  strength,  initial  elec¬ 
tron  density  profile,  etc.,  were  chosen  as  in  the  experiment  [69]. 
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Fig.  19  Initial  condition  for  wave  breaking  simulation,  (a)  Particle  curves,  (fo) 
Electron  density.  The  small  dark  gray  patch  is  four  times  as  dense  as  the  large 
background  patch. 
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Eig.  2U  Wave  breaking,  one  rotation,  (a)  Particle  curves,  (b)  Electron  density, 
Color  indicates  density  (red  is  high;  blue  is  low).  The  small  dark  gray  patch  is 
tour  times  as  dense  as  the  large  background  patch. 


A  small  dot  of  radius  is  0.1035  cm  and  high  density  nf\ni  = 
1  *0  x  11) J  cm”3  is  superimposed  on  a  large  background  electron 
disk  of  radius  0.8  cm  and  low  density  nt u*k  —  1*8  x  I0G  cm-3. 
Fig.  1 9  shows  the  initial  condition  with  (a)  curves  on  which  the 
particles  were  placed,  and  (b)  electron  density  mapped  onto  a 
mesh.  The  dot  rotates  about  the  center  ol  the  disk;  Fig,  20  shows 
the  solution  after  one  rotation  and  Fig,  21  show's  the  solution 
after  four  rotations.  As  the  del  rotates,  it  entrains  material  from 
the  disk  and  causes  a  wave  to  form  on  ihe  disk  boundary.  The 
results  are  in  excellent  agreement  with  experiment  [69],  Note 
lhal  the  fiiamentation  in  Fig.  21(a)  gives  rise  to  a  diffuse  region 
when  mapped  to  a  mesh,  as  in  Fig,  2 1(b),  Since  the  experimental 
results  are  obtained  by  crashing  the  electrons  onto  a  cathode 
ray  tube  and  capturing  the  image  with  a  charge-coupled  device 
camera,  the  resulting  spatial  blurring  may  be  similar  to  that  of 
mapping  the  particles  onto  a  mesh.  Hence,  the  diffuse  region 
seen  in  the  experiment  may  be  due  to  li lamentation  beyond  the 
experimental  resolution.  A  detailed  analysis  of  the  results  is  in 
preparation  [41]. 

IV.  Future  Work 

Our  short  term  goal  is  to  apply  BIT  and  particle  insertion  to 
simulate  the  warm  two-stream  instability.  In  addition,  we  are  ex¬ 
tending  our  2-D  field  solvers  to  three  dimensions,  as  well  as  op¬ 
timizing  the  ircecode  approach  for  solving  the  linear  system  as¬ 
sociated  with  ihe  homogeneous  solution.  In  the  future,  we  plan 
to  incorporate  a  mcsh-lree  DSMC  code  currently  under  develop¬ 
ment  to  permit  the  mesh-free  simulation  of  collisional  plasmas. 
We  plan  U>  begin  ihe  development  of  a  nonsiatistical  collision 
operator  which  will  accommodate  panicle  insertion.  An  addi- 
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Fig.  21  Wave  breaking,  four  rotations,  (a)  Particle  curves,  (b)  Electron 
density.  Color  indicates  density  (red  is  high;  blue  is  low)  The  small  dark  gray 
patch  is  four  times  as  dense  as  the  large  background  patch. 

lional  goal  is  to  incorporate  cl  lister-cluster  approximations  sim¬ 
ilar  in  spirit  to  the  fast  multipole  method  [27].  Finally,  we  plan 
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k)  investigate  whether  this  grid-free  approach  can  be  extended 
to  problems  with  lime- varying  magnetic  fields, 

V.  Conclusion 

We  have  applied  a  grid-free  boundary  integral/treecode 
(BIT)  field  solver  to  several  bounded  plasma  problems  and  ob¬ 
tained  comparable  or  belter  results  than  traditional  mesh-based 
methods.  We  demonstrated  that  BIT  exhibits  comparable 
timing  for  a  given  accuracy  and  is  capable  of  handling  com¬ 
plex  geometry,  as  well  as  a  mixture  of  boundary  conditions. 
Further  validation  is  underway.  We  believe  that  BIT  offers  an 
attractive  alternative  to  mesh-based  approaches  for  electrostatic 
problems.  The  extension  of  these  grid-free  methods  to  plasma 
problems  involving  time-varying  magnetic  fields  is  still  a 
significant  challenge. 
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